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We use determinant quantum Monte Carlo to study the single particle properties of quasiparticles 
and phonons in a variant of the two-dimensional Holstein model that includes an additional non¬ 
linear electron-phonon {e-ph ) interaction. We find that a small positive non-linear interaction 
reduces the effective coupling between the electrons and the lattice, suppresses charge-density wave 
(CDW) correlations, and hardens the effective phonon frequency. Conversely, a small negative non¬ 
linear interaction can enhance the e-ph coupling resulting in heavier quasiparticles, an increased 
tendency towards a CDW phase at all fillings, and a softened phonon frequency. An effective linear 
model with a renormalized interaction strength and phonon frequency can qualitatively capture this 
physics; however, the quantitative effects of the non-linearity on both the electronic and phononic 
degrees of freedom cannot be captured by such a model. These results are significant for typical non¬ 
linear coupling strengths found in real materials, indicating that non-linearity can have a significant 
influence on the physics of many e-ph coupled systems. 

PACS numbers: 


I. INTRODUCTION 

The electron-phonon (e-ph) interaction plays an im¬ 
portant role in many syst ems including conventional met- 
als and s upe rconductors,^^ organic semiconductors, 
fullereneSjlSEl and a large number of transition metal 
oxides.l^lti^In general, the e-ph interaction induces a local 
distortion of the lattice surrounding a carrier, resulting 
in quasiparticles dressed by phonon excitations known as 
polarons. In some cases these lattice distortions can be 
large and tightly bound to the electron, such that the 
quasiparticle has a very large effective mass m* and van¬ 
ishing quasiparticle residue Z^^ln this limit, the quasi¬ 
particle is typically referred to as a small polaron. In 
most models, the crossover to the small polaron regime 
occurs for A > 1, where A parameterizes the (dimension¬ 
less) strength of the e-ph interaction. When A < 1 the 
carriers are still dressed by the lattice, forming large po¬ 
larons where the lattice distortions are spread out over 
many lattice constants. 

Almost all of our knowledge about the effects of the 
e-ph interaction has been obtained from linear models. 
The derivation of these models is standard. First, the 
e-ph interaction is expanded in powers of the atomic dis¬ 
placement. This is followed by a truncation of the expan¬ 
sion to linear order under the assumption of small lattice 
displacements. These same models, however, often pre¬ 
dict the formation of small polarons or charge-density- 
wave (CDW) phases for sufficiently large e-ph coupling, 
which are characterized by large lattice displacements. 
For example, the displacements in the linear Holstein and 
Hubbard-Holstein models can be on the order of the lat¬ 
tice constant when CDW correlations are signihcant.l^^lt^ 
This can occur even for weak values of the e-ph coupling 
if the Fermi surface is well nested, as is the case for the 
Holstein model with nearest neighbor hopping on a cu¬ 
bic lattice. Similarly, a small polaron can be dressed 


by tens to hundreds of phonon quanta,^ implying the 
presence of a heavily distorted lattice surrounding the 
carrier. Clearly, these predictions directly violate the 
assumptions underlying the linear model, which is an 
unambiguous sign that important physics has been dis¬ 
carded during its derivation.!^ 

These considerations show that the higher order terms 
in the e-ph interaction are likely to be important when¬ 
ever the linear term is large (or when strong nesting con¬ 
ditions are present). But non-linear and anharmonic ef¬ 
fects can also be “switched on” in a weakly coupled sys¬ 
tem, if the underlying atoms of the lattice are driven far 
from their equilibrium positions by an external perturba¬ 
tion. For example, several recent experiments have ex¬ 
ploited optical pump pulses to drive the lattice, creating 
large lattice deformations or exciting coherent phonon 
oscillationsThe se ex cited lattice states can melt 
various ordered phase^^^ or can induce transient super¬ 
conducting states.!22EI] Here, nonlinear and anharmonic 
phonon dynamics are thou ght to p lay a vital role in cre¬ 
ating such transient states.^^SESEi] 

The first attempts to include higher order terms were 
made by Adolphs and Berdu (Ref. [211). They consid¬ 
ered the effects of non-linear e-ph interactions on a single 
polaron in the Holstein model using the momentum av¬ 
erage approximatioiP^ and found that small non-linear 
couplings dramatically undress the quasiparticle. This 
was attributed to a simultaneous hardening of the bare 
phonon frequency and a renormalization of the bare e- 
ph coupling constant, resulting in an overall weaker effec¬ 
tive linear interaction. Later work by some of the present 
authors considered finite carrier concentrations and tem¬ 
peratures using non-perturbative determinant quantum 
Monte Carlo (DQMC).!^Here too, the presence of a non¬ 
linear interaction was found to suppress the tendency to¬ 
wards the formation of CDW and superconducting states 
found in the linear model. 
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FIG. 1: (color online) The momentum dependence of the 
charge susceptibility xc(q) as a function of non-linear inter¬ 
action strength ^ > 0. The trends shown here follow those 
reported in Ref. [201 where the increasing non-linear interac¬ 
tion strength suppresses the CDW correlations in the system. 
The parameters for the calculation are as indicated. Error 
bars smaller than the marker size have been suppressed for 
clarity. 


It is important to note that the effects of the non-linear 
interactions observed in Refs. [50] and [53] are significant, 
even for relatively small non-linear interaction strengths. 
This is illustrated in Fig. which shows the suppres¬ 
sion of the CDW correlations in the half-filled Holstein 
model as a function of the non-linear coupling strength. 
[The degree of non-linearity is indicated by the ratio of 
the pre-factors of the quadratic (_g 2 ) and linear (gi) in- 

Here, the charge 


HA 


teraction terms £ = —, see Sec. 

^ . 9l’ 

susceptibility provides a measure ot the strength of the 
CDW correlations. It is defined as 


1 

Xc(q) = ^ dT{Trp{q,T)pH(i,0)), (I) 

where p(q) = ^ ^ and TV is the time-ordering 

operator. In this case, the well-knowiP^ES! tendency for 
the linear model to form a Cl = {TT/a,TT/a) CDW at half¬ 
filling is suppressed for £ as small as ~ 0.02 — 0.05. 

The suppression of xc(q) for small £ is noteworthy be¬ 
cause this ratio is typical of many models used to parame¬ 
terize e-ph interactions in real materials. In the transition 
metal oxides, for example, electrons can couple strongly 
to oxygen bond stretching modes via the modulation of 
the near-neighbor transition metal 3d-oxygen 2p (TM-0) 
hopping integ ral tp^, which depends on the TM-0 bond 
distance. b^F4T jf 

d is TM-0 bond distance and do is its 
equilibrium value, then this coupling mechanism leads to 
terms in the Hamiltonian of the form 

^kin = ^ ^ ^pdid)c\^^Pj^^, (2) 

(ij'l.o’ 


dependence tpd ~ {d/do)~^, with /? = 3.5.1^ Expanding 
this dependence to second order in powers of {d — dg) 
gives 

iTe-ph = ^ [o:i{d - do) + a 2 {d - do)^] (3) 

where ai = — f3tpd{do)/do and a 2 = ■ (The 

zeroth order terms are included in the non-interacting 
terms of the Hamiltonian). Taking typical value^^ for 
the various parameters in transition metal oxides, one 
arrives at a ratio of j^j ^ 0.05. Therefore real materi¬ 
als have intrinsic non-linear interactions that are large 
enough to be relevant when the linear coupling is strong. 
For this particular coupling mechanism the sign of £ is 
negative, however, this is not guaranteed for all coupling 
mechanisms.^ 

In this paper we expand upon our previous work exam¬ 
ining the impact of non-linear interactions on the CDW 
and superconducting phases of the Holstien model.l^ 
We present results for the single-particle electronic and 
phononic properties of the model, thus providing a more 
comprehensive picture of the effects of non-linearity. 
Since DQMC is formulated in the grand canonical en¬ 
semble, we are able to examine these properties at fi¬ 
nite temperatures and carrier concentrations for the first 
time. We find that the inclusion of a non-linear interac¬ 
tion renormalizes both the effective frequency of the Hol¬ 
stein phonon and the effective e-ph coupling strength, 
resulting in significant changes in both the electronic 
and phononic properties of the model. Furthermore, we 
demonstrate that while the qualitative effects of the non¬ 
linearity can be framed in terms of an effective linear 
model, the quantitative effects on both the electronic and 
phononic properties of the model cannot be. This conclu¬ 
sion requires an examination of both the electronic and 
phononic properties, and thus cannot be arrived at by 
considering electronic properties only. Our results reen¬ 
force the notion that the full non-linearity must be in¬ 
cluded in order to obtain an accurate picture of both 
the electronic and phononic degrees of freedom whenever 
strong linear e-ph interactions are present. 


II. METHODS 

A. The Non-linear Holstein Model 

We study a variant of the Holstein model that includes 
additional non-linear interaction terms. The Hamilto¬ 
nian is partitioned as 


iT = iTei + iTi,t + 7Ti,t, (4) 


iTel 




where the {Pj,<y) operators act on the TM 3d and O 2p 
orbitals, respectively. The hopping integral has a typical 


where 


2,(7 


( 5 ) 
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contains the non-interacting electronic terms, 
-ffiat = ^ 

i 

contains the non-interacting lattice terms, and 


p2 

% 


MVL 


-xi 


= E^ 


h\K + 2 


( 6 ) 


i?i„t = E = E (7) 

i,k,(T i,k,<7 


contains the interaction terms to order in the atomic 
displacement. Here, cj^ {Ci creates (annihilates) an 

electron of spin cr on lattice site i; b] (b-) creates (an¬ 
nihilates) a phonon on lattice site i; ^ 
is the number operator; /r is the chemical potential; t 
is the nearest-neighbor hopping integral; M is the ion 
mass; H is the phonon frequency; Xi and Pi are the lat¬ 
tice position and momentum operators, respectively; and 
gk = ak{ 2 MQ )~2 is the strength of the e-ph coupling to 
order in displacement. 

The non-linear Holstein model is characterized by sev¬ 
eral dimensionless parameters, and the specific choice 
in parameterization is not unique. Her e, we follow the 
convention used in previous works ,^22123] the usual 

dimensionless parameter A = a\/{MVPW) = g\/{4i^) 
parameterizes the linear coupling strength and = 
gk/gk-i parameterizes the non-linear interaction terms. 
This choice provides a convenient interpretation with 
large A implying a strong linear interaction and large 
implying strong non-linear effects. In the linear model 
(Cfc = 0) A > 1 implies the formation of small po- 
larons. Thus this choice of parameterization is also useful 
for making comparisons to our expectations gained from 
studying the linear model. 


B. Determinant Quantum Monte Carlo 


therefore able to perform simulations to arbitrarily low 
temperatures,^ however, we find that most of the physi¬ 
cal properties we are interested in here can be examined 
for j3 = A/1. We use this temperature for all plots in this 
work unless stated otherwise and present results for an 
imaginary time discretization of At = 0.1/t. In all of 
our simulations we have not observed any significant At 
errors introduced by this choice. 


C. Analytic Continuation 

The DQMC calculation provides the phonon Green’s 
function Z3(q, t) = (TXq(T)A_q(0)) measured on the 
imaginary time axis. The phonon Green’s function is 
related to the phonon spectral function on the real axis 
by the integral equation 


f dTD{ci,T)= r 
Jo J-i 


dcu B{c[,oj) 
2tt uj 


( 8 ) 


This equation also provides a normalization condition for 


H(q, w). In section HIG we examine the phonon spec¬ 


tral properties of the non-linear Holstein model, which 
requires that the phonon Green’s function be analyti¬ 
cally continued to the real frequency axis. This is accom¬ 
plished with the Maximum Entropy method (MEM).Sil 
The analytic continuation procedure is identical to the 
one given in Ref. HU and the reader can refer to there for 
details. 

MEM requires a model default function to define the 
entropic prior. We adopt a momentum-independent 
Lorentzian model, which is peaked at the renormalized 
phonon frequency predicted by the mean-held treatment 
of the non-linear interaction VLmf = H -|- 2 g 2 (see sec¬ 
tion IhT^ and of width t. Our results are robust against 


reasonable changes in this choice of default function. 


We use DQMG to study the non-linear Holstein model. 
The details of the method are given in several references 
(see for example Refs. [371-[M]) and the specihcs for han¬ 
dling the lattice degrees of freedom can be found in Refs. 

muni andm 

In our calculations we keep gi > 0 without loss of 
generality. Furthermore, Ref. |33| examined terms to 4*** 
order in the interaction and found that the largest ef¬ 
fect was produced by the 2'^'^ order terms. We expect a 
similar result here and restrict ourselves to k = 2 while 
dehning C = If"- Throughout this work we examine two- 
dimensional square lattices with a linear dimension N (a 
total oi N X N sites) and set a = t = M = la.s our units 
of distance, energy, and mass, respectively. We typically 
work on lattice sizes ranging from iV = 4 to 8 in size. In 
general we do not observe signihcant hnite size effects,!^ 
which is likely due to the local nature of the interaction 
in the model. 

The Holstein model and its non-linear exte nsion do 
not suffer from a fermion sign problem !^^ l ^°l We are 


III. RESULTS AND DISCUSSION 

Our results for the single-particle electronic and 
phononic properties of the non-linear Holstein model are 
presented in this section. We will begin by first exam¬ 
ining the renormalization of the quasiparticles as a func¬ 
tion of doping, temperature, and phonon energy. We 
then examine how the interplay between the quasiparti¬ 
cles and renormalized phonons affect the energetics of 
the system. Following this, results are presented for 
the renormalization of the phonons. Finally, after all of 
these effects are examined, we demonstrate that the si¬ 
multaneous renormalization of the electronic and phonic 
subsystems cannot be quantitatively captured by an ef¬ 
fective linear model. Gombined, these results paint a 
more detailed picture of how both the quasiparticles and 
phonons are renormalized by the non-linear e-ph interac¬ 
tion, which cannot be obtained by examining only one of 
these subsystems. In the following sections we consider 













4 



FIG. 2: (color online) The quasiparticle residue ^(k) as a 
function of band filling (n) for (a) k = (0, 0), (b) (0,7r/2), (c) 
(0,7r), (d) (7r/2,7r2), (e) (7r/2,7r), and (f) (tTjTt). Results are 
shown for various values of the non-linear interaction strength 
as indicated in panel (f), and are obtained using an Ai = 
4x4 cluster with a linear coupling A = 0.25 and an inverse 
temperature /3 = 4/t. Error bars smaller than the marker size 
have been suppressed for clarity. 


results for the case ^ > 0. Considerations of the ^ < 0 
case are presented in Sec. [ml 


A. The quasiparticle residue 


We begin by examining the polaron’s quasiparticle 
reside ^(k) as a function of the non-linear coupling 
strength and doping. This quantity is related to the ef¬ 
fective mass via Z~^ cx It can be obtained from the 
imaginary axis self-energy E(k, using the relation¬ 
ship Z(k) = where 


6(k) 


lim 


5S'(k,fa;„) 

duJn 


( 9 ) 


Here, we approximate 6(k) by evaluating Eq. (j^ for the 
lowest Matsubara frequency a;„ = n 1/3. 

Fig. i shows Z{\i) as a function of carrier concen¬ 
tration for several values of the non-linear coupling 
These results were obtained on an TV = 4 cluster, using 
a linear coupling strength A = 0.25 and il = t. In the 
linear model = 0, red dots) the quasiparticle residue 
decreases as the filling approaches (n) = 1, where the 
Q = (tt, tt) CDW correlations begin to dominate the sys¬ 
tem (Fig. [^. Note that strong CDW correlations are 
observed, even for the small value of the linear coupling 



FIG. 3: The spectral weight at the Fermi level given by 
l3G{r = 0,T = (5/2) = (3G^/2 as a function of band filling 
in) for various values of the non-linear coupling strength C 
(a) Results for a A = 4 cluster and an inverse temperature 
d = 4/t. (b) Results for a larger N = 8 cluster and a lower 
temperature (J = 5/t. All results are obtained for a linear 
coupling A = 0.25 and a frequency of = f. Error bars 
smaller than the markers have been suppressed for clarity. 


used here, due to a perfect (tt, tt) nesting condition in 
the two-dimensional Fermi surface. This nesting condi¬ 
tion also results in strong lattice displacements in the 
linear model. As a result, the inclusion of the non-linear 
terms has a significant effect on the quasiparticle residue 
where, for ^ > 0, a significant undressing of the polarons 
occurs and the quasiparticle residues at all momenta be¬ 
gin to rise. This occurs at all doping, however, the effect 
is more pronounced near half-filling. (Our C = 0 results 
are in good agreement with Ref. 1441 which examined 
larger system sizes using a complementary diagrammatic 
Monte Carlo method.) 

The formation and suppression of the CDW gap is also 
reflected in the spectral weight at the Fermi level, which 
can be obtained from the local imaginary time Green’s 
function /3G(r = 0,r = [3/2) = /3G,g/2- Fig. plots 
G/ 3/2 as a function of filling (n) for the same parameters 
used in Fig. H Fig. ir plots similar results obtained 
on a larger cluster and at lower temperature, where the 
qualitative behavior is the same. The spectral weight 
in the linear model initially grows with increasing car¬ 
rier concentration, but saturates as the concentration 
approaches half-filling and CDW correlations begin to 
dominate. When a non-linear interaction is introduced, 
however, G ^/2 increases at most fillings, which is most 
pronounced near (n) ^ 1. (The dip around (n) = 0.6 
is a finite size effect due to the smaller number of mo¬ 
mentum points in the N = 4 cluster. It is much less 
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FIG. 4: (color online) A finite size scaling analysis of Z(0, tt) 
as a function of 1/A^ where N is the linear dimension of the 
cluster. The parameters for the calculations are {3 = 4/t, 
Q = t, and A = 0.25. Error bars smaller than the marker size 
have been suppressed for clarity. 


pronounced on the larger N = 8 cluster.) This spectral 
weight increase directly reflects the increase in the quasi¬ 
particle residue and the suppression of the CDW correla¬ 
tions. Previously we showed that a large non-linear cou¬ 
pling drives the system into a metallic state at half-filling, 
with the value of (3Gp/2 approaching the non-interacting 
value.l^ The results in Fig. indicate that this also oc¬ 
curs for carrier concentrations away from half-filling. 

The results presented in Figs. and are obtained 
on a iV = 4 cluster; however, they are qualitatively repre¬ 
sentative of the results obtained for all examined cluster 
sizes, as hinted at by comparing Figs. an d[^. To con¬ 
firm this, in Fig. [^we perform a finite size scaling analy¬ 
sis for Z(0, tt) at half-filling, where the reduction in Z by 
CDW correlations is most pronounced. From this analy¬ 
sis it is clear that the qualitative behavior is not affected 
by finite size effects and survives in the thermodynamic 
limit. Moreover, the more pronounced finite size effects 
occur when the non-linear interaction is weak. As similar 
scaling results were obtained for both the charge suscep¬ 
tibility and the electron spectral weight in our previous 
work,l2Sl we conclude that the qualitative physics of the 
non-linear model can be obtained on an = 4 cluster. 

In Fig. we consider the temperature dependence 
of 2'(k) and average displacement of the lattice {X) = 
J2i=i ^i,i half-filled model. Here, re¬ 

sults for Z(0, tt) only are shown, since similar trends were 
found at all momenta. Focusing first on the linear model, 
we find that Z(0,7r) decreases with temperature as the 
CDW correlations begin to set in. The average lattice 
displacement, however, does not exhibit the same tem¬ 
perature dependence (see inset of Fig. [^).l^ As the non¬ 
linear interaction strength grows, however, the quasipar¬ 
ticle residue increases back towards its non-interacting 
value. For small ^ this rise is somewhat rapid, but it gives 
way to a more gradual increase for ^ ^ O.I. The increase 
is also accompanied by a decrease in the average lattice 
displacement (inset of Fig. [^. This behavior mirrors 
the observed ^-dependence of the CDW susceptibility,!^ 
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FIG. 5: The (a) temperature and (b) Q dependence of the 
quasiparticle residue in the half-filled model as a function of 
non-linear coupling strength C The insets show the corre¬ 
sponding expectation value of the lattice displacement. All 
results are obtained on an A = 4 cluster and with a linear 
coupling A = 0.25. Error bars smaller than the marker size 
have been suppressed for clarity. 


and is consistent with the conclusion that a finite ^ > 0 
undresses the carriers and relaxes the lattice distortions 
normally present in the linear model. 

The ri-dependence of Z{tt,0) and {X) for the same 
model are shown in Fig. Here, the ^ = 0 results 
are consistent with those obtained for the ID Holstein 
model, where the tendency to form a CDW grows with 
decreasing phonon frequencies.!^ Consequently, both the 
quasiparticle residue and average lattice displacement de¬ 
crease as the value of increases. The introduction of 
^ > 0 results in the further decrease in these quantities. 

From this section we conclude that the non-linear in¬ 
teraction with ^ > 0 acts to undress the quasiparticles 
and that this is a generic result, regardless of the values 
of and /3. The undressing, however, is much more pro¬ 
nounced at low temperatures, for smaller values of the 
phonon frequency, and near half-filling, where the CDW 
correlations (and subsequently the local lattice displace¬ 
ments) are largest in the linear model. 


B. Electron and Phonon energetics 

The average kinetic energy of the electronic subsystem 
(KE)ei = at /3 = 4/t is shown in Fig. 

1^ Results are shown as a function of band filling (n) for 
various ^ and for a linear coupling A = 0.25. For ^ = 0 
the total kinetic energy — (KE)e increases as a function 
of (n) as higher momentum states are populated in the 
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FIG. 6: (color online) (a) The electron kinetic energy (KE)ei 
and (b) the e-ph energy (E)e-ph as a function of band hlling 
(n) for various non-linear interaction strengths as indicated. 
Results are obtained on an = 4 cluster and with a linear 
coupling A = 0.25, phonon frequency = t, and an inverse 
temperature (3 = 4/t. The dashed lines in panels (a) and (c) 
indicate the non-interacting result, (c) and (d) show corre¬ 
sponding results for a larger N = 8 cluster and 13 = 5/t. The 
remaining parameters are the same as in (a) and (b). Error 
bars smaller than the marker size have been suppressed for 
clarity. 


Fermi sea, however, the total kinetic energy saturates 
as the filling increases beyond (n) > 0.6. This is due 
to the saddle point in the band dispersion at (0,7r) and 
is also present in the non-interacting model (indicated 
by the dashed line). When the non-linear interaction 
is added we see an overall increase in the total kinetic 
energy, which tends towards the non-interacting value at 
all fillings for large This again reflects the undressing of 
the quasiparticles and the subsequent increase in mobility 
of the electronic subsystem. 

Fig. [§D shows the corresponding e-ph interaction en¬ 
ergy, defined as (E)e-ph = g 2 htXf). Unsur¬ 

prisingly, the total e-ph interaction energy increases with 
band filling as both the average number of electrons per 
site and the average lattice displacement increase. This is 
most evident in the linear model (^ = 0). Increasing the 
value of ^ naturally leads to smaller lattice displacements 
and a significant decrease in (E)e-ph- 

The average kinetic (KE)ph and potential (PE)ph en¬ 
ergies of the lattice are shown in Eigs. [3^ and[^, respec¬ 
tively. They are given by 


(KE)ph 


(PE)ph 



( 10 ) 

( 11 ) 


(The factor of appearing in Eq. (10) is a Eu¬ 
clidean correction introduced by the Wick rotation to 


FIG. 7: (color online) The phonon (a) kinetic energy (KE)p(, 
and (b) potential energy {PE}ph as a function of band filling 
(n) for various of the non-linear interaction strength as 
indicated in panel (b). Results are obtained on an A'’ = 4 
cluster and with a linear coupling A = 0.25, = l,and f3 = 

4/t. (c) and (d) show similar results obtained on a larger 
N = 8 cluster with j3 = 5/t. Error bars smaller than the 
marker size have been suppressed for clarity. 


the imaginary-time axis.^1^ 

In the linear model we see a very weak variation in 
the phonon kinetic energy as a function of filling, with a 
slight decrease observed near half-filling when the CDW 
correlations increase. This is consistent with prior obser¬ 
vations of the lattice kinetic energy in the viciniW of a 
CDW transition in the Hubbard- Holstein model.^ The 
average potential energy of the lattice grows as the av¬ 
erage number of carriers per site increases. When the 
non-linearity is introduced and the lattice distortions di¬ 
minish (see inset), we see an increase in the lattice 
kinetic energy, which is attributed to the hardening of the 
phonon dispersion. At the same time, we see a decrease 
in the total lattice potential energy. Here, the non-linear 
interaction has two opposing effects: the increase in the 
phonon frequency increases the lattice potential energy 
while the decrease in the effective linear coupling decrease 
the net lattice distortions and subsequently lowers the po¬ 
tential energy. Our results indicate that the latter effect 
has the stronger impact. 

The energetics reported here are completely consistent 
with the conclusion that the non-linear interaction acts 
to harden the phonon frequency and weaken the effective 
linear interaction, which results in an undressing of the 
Holstein polaron for ^ > 0. 


C. Phonon Spectral Properties 


In the linear Holstein model the formation of the CDW 
phase is accompanied by a softening of the phonon dis¬ 
persion to zero energy at the nesting wavevector Q = 
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FIG. 8: (color online) The phonon density of states Afph(i') = 
1 /) for the half-filled model as a function of the 
non-linear interaction strength. The remaining calculation 
parameters are as indicated. 


2'kp = (tt, 7r).l22lllHini This softening is generated by the 
strong nesting condition of the non-interacting Fermi sur¬ 
face. The inclusion of the non-linear e-ph interaction is 
therefore expected to modify the phonon dispersion in 
two important ways: first, it will undo the softening at 
Q as the CDW correlations are suppressed. Second, it 
will result in an overall renormalization of the phonon 
frequency, as discussed in the introduction. We con¬ 
firm these expectations in this section by examining the 
phonon spectral function i?(q, v) and the phonon density 
of states (DOS) Nph{v) = ^ ^)- 

In the non-interacting limit, i3(q, z^) and Nph(y) are 
delta functions centered at the bare phonon frequency 
O. In the presence of a non-zero linear interaction only, 
this distribution shifts to lower energy and broadens. 
This is illustrated in Figs. |^and 1^, which plot Nph{i/) 
and the momentum-resolv^ phonon spectral function 
respectively, for the half-hlled model. These 
results were obtained on IV = 8 clusters, with a linear 
coupling A = 0.25, = t, and at an inverse temperature 

of /3 = 4/t. Due to the finite value of A, the phonon fre¬ 
quency softens from its non-interacting value and Nph{v) 
for the linear model consists of a broad, asymmetric dis¬ 
tribution centered at ~ 0.60t. The asymmetry in Nph{i') 
reflects the momentum dependence of the softening and 
the low-energy spectral weight in B(Q,i/), coupled with 
the requirement that B{0) — 0 for bosons. This is more 
easily seen in the momentum-resolved spectral function 
(Fig. ^), which has a clear Kohn anomaly at Q = (tt, tt). 

Two prominent changes occur when ^ 7^ 0. First, the 
peak in the DOS shifts to higher energies, which verifies 
the hardening of the effective phonon frequency. This 
behavior is also clearly seen in the momentum resolved 
spectral functions, shown in Fig. Second, the pro¬ 
nounced Kohn anomaly begins to disappear as the CDW 
correlations are suppressed with increasing values of 
Both of these results confirm our expectations. 
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FIG. 9: (color online) The momentum resolved phonon spec¬ 
tral function B{q, v) for the half-filled model and for various 
values of the non-linear interaction strength, as indicated in 
each panel. Results were obtained on an = 8 cluster with 
A = 0.25, D = t, /3 = 4/t and Ar = 0.1/t. The black squares 
indicate the position of the peak in the phonon spectral func¬ 
tion. 



D. Mean-field Treatment of the quadratic 
e-ph interaction 

As we have repeatedly seen, the non-linear e-ph inter¬ 
action acts to renormalize both the bare linear interac¬ 
tion strength A and the bare phonon frequency Q. Both 
of these effects can be qualitatively understood at the 
mean-held (MF) level for the quadratic model, where an 
effective linear Hamiltonian is obtained by performing a 
MF decoupling of the interaction terms proportional to 
h\b\ and bibi^ The resulting effective MF Hamiltonian 
is 


Hmf = Bel + ^MF (blb^ ^ j 

+ ^9MFn^,a { 4 -hb-'j , ( 12 ) 

i,<T 

where ^mf = D -f 2^2 and qmf = 5i(l - are 

the renormalized phonon frequency and e-ph coupling 
constants, respectively. One immediately sees that the 
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FIG. 10: (color online) A comparison of (a) the quasiparti¬ 
cle residue and (b) CDW and pair field susceptibili¬ 

ties xsc obtained from the non-linear model and its effective 
linear model, as defined in the main text. The bare linear 
coupling and phonon frequency are A = 0.25 and Q. = t, re¬ 
spectively. In both cases results are obtained on an = 8 
cluster and at an inverse temperature of /3 = 4/t. Error bars 
smaller than the marker size have been suppressed for clarity. 


quadratic e-ph interaction leads to a softening (harden¬ 
ing) of the phonon frequency and an increase (decrease) 
in the effective linear interaction strength gi for ^ < 0 
(^ > 0). These two effects combine to produce an overall 
increase (decrease) in the strength of the effective dimen¬ 


sionless coupling Ae// oc 

The MF treatment of the non-linear interaction is con¬ 
sistent with the general trends reported here and in Refs. 
[20land[2^ We stress, however, that the MF description 
only provides a qualitative picture of the non-linear ef¬ 
fects. To illustrate this, we compare our DQMC results 
for the full non-linear Hamiltonian against the predic¬ 
tions obtained from two sets of effective linear models. 
The first is the MF-derived model defined by Eq. (12). 
The second is the set of effective linear models whose pa¬ 
rameters are obtained by tuning the fleff and geff to 
reproduce the electronic properties of the system. 

We consider the MF-derived model first. Fig. [T^ com¬ 
pares the results for the quasiparticle residue, Xc(’’'i'^)j 
and the pair-field susceptibility xsc calculated using the 
full non-linear model [Eq. @] to results obtained from a 
DQMC simulation of the corresponding MF-derived lin¬ 
ear model [Eq. @] at half-filling. We find that the MF 
model does a poor job in quantitatively capturing the 
electronic properties; it underestimates both the quasi¬ 
particle residue Z{0, tt) and the tendency towards the for¬ 
mation of a CDW when compared to the full non-linear 
model. The MF model also over-predicts the magnitude 



FIG. 11: (color online) A comparison of the results obtained 
for the full non-linear model and an effective linear model 
where the value of the e-ph coupling constant has been ad¬ 
justed to reproduce the electronic properties of the non-linear 
model. Panels (a)-(d) show a comparison for an effective lin¬ 
ear model with D = t, equal to the bare phonon frequency. 
Panels (e)-(f) show a comparison for an effective linear model 
with fl = Hmf- The top row [panels (a) & (e)] compares the 
quasiparticle residues obtained with both models. The sec¬ 
ond row [panels (b) & (f)) shows the resulting charge (solid 
lines) and pair-field (dashed lines) susceptibilities. The third 
row [panels (c) & (g)] show the resulting phonon potential 
and kinetic energies. The potential energy has been divided 
by a factor of three and is indicated by the solid lines while 
the kinetic energy is indicated by the dashed lines. Finally, 
the bottom row [panels (d) & (h)) show the average value of 
the lattice displacement. The remaining parameters of the 
simulation are /3 = 4/t and N = 8. 


of the superconducting pair-field susceptibility when ^ is 
large and under-predicts it when ^ is small. 

The results shown in Fig. demonstrate that the 
MF treatment of the quadratic interaction can only pro¬ 
vide a qualitative picture of the physics of the non-linear 
model; however, another choice in effective model might 
do a better job. We explored this possibility by adjusting 
the effective coupling strength in the linear model such 
that the linear model reproduce the electronic properties 
of the full non-linear model. This procedure was per¬ 
formed for two choices in the phonon frequency. First, 
we set the phonon frequency equal to the bare value and 
adjusted the value of the coupling strength to reproduce 






















































9 
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Momentum [1/a] 

FIG. 13: (color online) The momentum dependence of the 
charge susceptibility xIq) S'® a function of non-linear interac¬ 
tion strength ^ < 0 at half filling. The parameters are set 
as: A = 0.25, = 2t, /3 = 4. The results are obtained on an 
N = 8 cluster. Error bars smaller than the marker size have 
been suppressed for clarity. 


FIG. 12: (color online) (a) the quasiparticle residue and (b) 
CDW xc(7r, tt) (open symbols) and pair field susceptibilities 
Xsc (solid symbols) as a function of band filling. The param¬ 
eters are set as: A = 0.25, 12 = 2t, /3 = 4. The results are 
obtained on an = 8 cluster. Error bars smaller than the 
marker size have been suppressed for clarity. 


the phononic and electronic properties of the system. A 
similar conclusion was reached in Ref. |13] in the single 
carrier limit. 


E. Negative values of ^ 


the quasiparticle residue, as shown Fig. [n}i. The value 
of the linear coupling strength ()e// needed to produce 
this agreement is shown in the inset (black solid A), 
where it is compared against the corresponding value of 
gMF = gi + 2(72- By tuning the value of geff we are 
able to accurately capture the quasiparticle residue. The 
charge and superconducting pair-field susceptibilities are 
also well reproduced, indicating that this effective model 
is capable of capturing the electronic properties of the 
system. But when we examine the phonon properties 
(Figs. 11: & we find some disagreement, partic¬ 
ularity with respect to the predicted kinetic energy of 
the lattice, where the linear model systematically under¬ 
predicts the correct results. 


The comparison between the two models can be im¬ 
proved somewhat if we set the phonon frequency to be 
equal to Hmf and again readjust the value of geff - This 
case is shown in Fig. [15 -h. Using this choice we are 
again able to accurately capture the electronic properties 
and improve the comparison between the kinetic energy. 
But this comes at the expense of the level of agreement 
between the average lattice potential energy and the av¬ 
erage lattice displacement. From this we conclude that 
an effective linear model cannot capture both the elec¬ 
tronic and phononic properties of the non-linear model, 
for a fixed value of the phonon frequency. These results 
show that while the qualitative effects can be understood 
using an effective linear model, the full non-linear in¬ 
teraction should be retained if one wishes to accurately 
capture the effects of the non-linear interaction on both 


We have shown that a positive (^ > 0) non-linear cou¬ 
pling results in a hardening of the phonon frequency and 
a renormalization of the effective linear e-ph coupling to 
weaker values. But what about the case when ^ < 0, 
where the MF model predicts an enhanced effective lin¬ 
ear coupling? Before examining this case, we note that 
a large negative ^ necessarily requires the inclusion of 
additional anharmonic terms in the lattice potential.!^ 
For ^ < 0 ((72 < 0) the phonon frequency given by 
rie// = it + 2 g 2 can become negative for sufficiently large 
values of 772 , indicating an instability in the lattice. In 
this event the anharmonic terms of the lattice potential 
are required to maintain stability. At present, our codes 
do not contain such terms and we are unable to examine 
this case in great detail. We therefore restrict ourselves 
to a larger value oi il = 2t and small values of \g 2 \ in 
order to get a feel for the (72 < 0 regime while ensuring 
the stability of the lattice. 

Fig. 12 1 shows the quasiparticle residue, xci^^, and 
Xsc as a function of band filling for various values of 
^ < 0. These results were obtained for a linear coupling 
of A = 0.25 and on an A = 8 cluster. We find that the 
quasiparticles are more effectively dressed when ^ < 0 , 
and the quasiparticle residue is much smaller for all fill¬ 
ings when increasing negative quadratic interactions are 
included. The CDW correlations are also significantly en¬ 
hanced, as reflected in the charge susceptibility shown in 
Fig. Both of these observations are in line with the 
expected increase in the effective linear coupling. Fur¬ 
thermore, since the CDW phase directly competes with 
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s-wave superconductivity, the pair-field susceptibility is 
suppressed at filling values where the CDW correlations 
dominate. In addition, we also see a noticeable decrease 
in the pair-field susceptibility at band fillings where the 
CDW does not dominate. This suggests that the negative 
non-linear interaction has a detrimental effect on the su¬ 
perconducting transition temperature, which stems from 
the decrease in the quasiparticle residue. Fig. [T^ plots 
the momentum dependence of the charge susceptibility 
x(q) as a function of negative ^ at half filling, where it is 
clear that the dominant CDW correlations are still being 
set by the Fermi surface nesting condition. 

IV. SUMMARY AND CONCLUSIONS 

We have examined the role of non-linear e-ph in¬ 
teractions in shaping the single-particle electronic and 
phononic properties of Holstein polarons at finite car¬ 
rier concentrations and temperatures using DQMC. We 
find that the inclusion of a positive non-linear interac¬ 
tion term serves to undress the polaron leading to carri¬ 
ers with lighter effective masses. This leads to changes in 
the energetics of both the electrons and phonons, as well 
as the relaxation of the local lattice distortions surround¬ 
ing each carrier. This is due to a simultaneous hardening 
of the phonon frequency and renormalization of the ef¬ 
fective linear coupling to smaller values. We have also 
examined the case when the quadratic e-ph interaction 
has the opposite sign as the linear interaction, although 
this case cannot be explored in detail without the inclu¬ 


sion of additional anharmonic terms in the lattice po¬ 
tential. Nevertheless, in our limited range of accessible 
parameters, we find that a quadratic interaction results 
in an increased dressing of the carriers and an enhanced 
tendency towards the formation of a Q = (tt, tt) CDW 
ordered phase. 

While many of the effects we have discussed can be 
understood qualitatively at the mean-field level, we have 
demonstrated that the quantitative effects can only be 
captured by the full non-linear model. Specifically, the 
effective linear models fail to simultaneously capture the 
electronic and phononic properties. Therefore the full 
non-linear model must be retained if one wishes to ac¬ 
curately capture the properties of the electrons and the 
phonons. Our results are in good agreement with the 
results obtained in the single particle limit ,1^ and show 
that non-linearities are relevant at finite carrier concen¬ 
trations. 
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